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ORIGINAL PAGE 



Amiram Harten 

— 1946 - 1994 


Amiram Harten, a professor of Mathematics at Tel-Aviv University and a consultant at ICASE, 
died of a massive heart attack on August 5, 1994. He was 47. He was associated with ICASE 
since 1976 and had a great impact on the ICASE program of research in numerical analysis and 
algorithm development. He was an active participant in ICASE activities and visited ICASE at 
least once a year. His last visit to ICASE was during the week of May 22, 1994, when he participated 
in the Parallel Numerical Algorithm Workshop. He planned to come back on August 29, 1994 
to be at ICASE for a month. His untimely passing away is an irreparable loss to ICASE and the 
mathematics community. He will be greatly missed by his friends and colleagues even as 

his influence lives on. 

This report is the written version of his lecture at the Parallel Numerical Algorithm Workshop, 

and it is sadly noted to be his last ICASE report. 
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Multiresolution Representation and Numerical Algorithms: 

A Brief Review 

Ami Harten 1 

School of Mathematical Sciences 
Tel-Aviv University 
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Department of Mathematics 
UCLA 

ABSTRACT 

In this paper we review recent developments in techniques to represent data in terms of its 
local scale components. These techniques enable us to obtain data compression by eliminating 
scale-coefficients which are sufficiently small. This capability for data compression can be 
used to i educe the cost of many numerical solution algorithms by either applying it to the 
numerical solution operator in order to get an approximate sparse representation, or by 

applying it to the numerical solution itself in order to reduce the number of quantities that 
need to be computed. 
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1. Introduction 


Fourier analysis, which provides a way to represent square-integrable functions in terms 
of their sinusoidal scale-components, has contributed greatly to all fields of science. The 
main drawback of Fourier analysis is in its globality - a single irregularity in the function 
dominates the behavior of the scale-coefficients and prevents us from getting immediate 
information about the behavior of the function elsewhere. 

The recent development of the theory of wavelets (see [Me] and [Ma]) was a great step 
towards local scale decomposition, and has already had great impact on several fields of 
science. In numerical analysis representation by compactly supported wavelets (see [Da] and 
[CDF]) is used to reduce the cost of many numerical solution algorithms by either applying it 
to the numerical solution operator to obtain an approximate sparse form (see [BCR]), or by 
applying it to the numerical solution itself to obtain an approximate reduced representation 
in order to solve for less quantities (see e.g. [LT], [MR] and [BMP]). The main drawback of 
the theory of wavelets is that it attempts to decompose any square integrable function into 
scale-components which are translates and dilates of a single function. Consequently there 
are conceptual difficulties in extending wavelets to bounded domains and general geometries. 
Furthermore, the uniformity of the underlying wavelet approximation makes it impossible to 
obtain an adaptive (data-dependent) multiresolution representation which fits the approx- 
imation to the local nature of the data. The only adaptivity which is possible within the 
theory of wavelets is through redundant “dictionaries.” 

In a series of works [HI -3] we have studied the question of how to represent discrete data 
which originates from unstructured grids in bounded domains in terms of scale decompo- 
sition. Combining ideas from multigrid methods, numerical solution of conservation laws, 
hierarchial bases of finite element spaces, subdivision schemes of Computer-Aided Design 
and of course - the theory of wavelets, we came up with the more general concept of “nested 
sequence of discretization.” Given discrete data which can be associated with a nested 
sequence of discretization we show that it has a multiresolution representation, i.e., a one- 
to-one correspondence between the given data and its scale-decomposition. This framework 
is a generalization of the theory of wavelets in the sense that under conditions of uniformity 
its natural result is wavelets. 

In this paper we review the work in [HY], [AC], [ACD], [H4-5] where the previously men- 
tioned works on numerical solution algorithms with representation by wavelets are extended 
to the more general framework of nested discretization. 
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2. Nested Discretization 


In this section we describe the class of discrete data for which we can obtain representation 
in terms of a scale-decomposition. We start with two examples. 

Example 1. 

Let us consider continuous functions / in the interval [0,1] 

fGf = C°[ 0,1], 


and let {X k }% =0 


be the following nested sequence of uniform grids 

X* = {^}So, 4 = ihk, h k = 2 J k = 2 k J 0 


(2.1a) 


for some integer J 0 with h 0 = 1/Jo- Here k = L is the finest grid and k = 0 is the coarsest. 
Observe that X k is obtained by dyadic refinement of X fc_I , i.e. 

4 = 44 i = 0,...,J fc -i (2.16) 

= i = 1 ,- Jfc-i- (2-lc) 

Let us consider now the discrete values v k = {4}/=o which are obtained from point value 
discretization of / € T in the k - th grid 

4 =: (V k f)i = f(x k ), i = 0, . . . , J k , 0 < k < L. (2.2 a) 

It follows from (2.1b) that v k ~ 1 can be obtained from v k by the decimation 

4~* = 4, i = 0,...,J fc _i. (2.26) 


Let I k ! (;r;i>* *) denote any continuous function in [0,1] which interpolates v k 1 at the 
grid points of X k ~ x , i.e. 

I k - 1 (x-v k - 1 )eJ r =C°{ 0,1] (2.3a) 

I k -^ x k-i. v k-i ) = v k-i for i = o (2.36) 

e.g. we can take 7 fc_1 (x; u* -1 ) to be the piecewise-linear interpolation 


I k - 1 (x-,v k ~' 1 ) = wf_" 1 1 + 


1 - v k :!), for *{:' < x < 4- 1 . 

hk - 1 


Given v k 1 we can approximate v k by 

4 » /*-*(*?; v*- 1 ), * = 0, . . . , Jfc. 


(2.4) 
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Let us denote 


e " = v--I k - 1 (x l ‘;v k - 1 ), t = 0, . . . , J fc (2.5a) 

and refer to e = as the prediction error. We observe that 

4 = 0 for i = 0,...,J k _ 1 (2.56) 

and denote the interpolation error at the odd grid points of X k by d k = {d k } J ] i~ l 1 

d 3 =: 4-1 = 4-1 ~ /fc_1 (4-i; j = • -,Jk- 1 . (2.5c) 

Clearly there is a one-to-one correspondence between v k and {d fc ,n* -1 }: Knowing v k we get 
v k ~ l and d k by (2.2b) and (2.5c), respectively. From {d k ,v k ~ 1 } we recover v k by 


V < 


2 i ~ v i 1 ) * — 0, ... , Jfc-1 

2i-l = I* 1 { x 2i- V V>C 1 ) + 4 1 * = 1, • • ■ , Jk-1- 


(2.6) 


one 


Since v k 1 can be likewise represented by {d k ~\v k ~ 2 } we get that there is a one-to- 
correspondence between the values of the finest level v L and the sequence of ( J L + 1) elements 
{d L i ■ ■ ■ ,d 1 , v 0 } which we denote by v 


L 1:1 


V <- 


+ {d L , . . . , rf 1 , v 0 } =: v M - 


(2.7) 

We refer to as the multiresolution representation of v L . It follows from (2.2b) and (2.5c) 
that the direct multiresolution representation (MR) transform v L v M can be expressed 
by the following algorithm 


f DO k = 


4 1 = 4, * = o,...,j fc _! 


(2.8a) 


d J = 4-1 - Ik ~ 1 (^ j -vV k ~ 1 ), j = 1,..., Jfc.,. 

We denote this transform by M , i.e. 


v M = M- v L . (2.8 b) 

Similarly it follows from (2.6) that the inverse MR transform v M <-► v L can be expressed by 


f DO k = 


4 = v-\ * = • • • j Jk—l 


(2.9a) 


1 4-i = /*- , (4-i;”*- 1 ) + 4 > = 
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and we denote 


and we denote 

v L = M -1 • VM- 

(2.96) 

We refer to d k as the scale-coefficients of the fc-th level of resolution, 
linear interpolation (2.4) we get from (2.5c) that 

For the piecewise- 

< = 4-i - ^‘-i 1 + 4 

(2.10a) 

in terms of / € T for which v k = V k f this can be expressed by 


d)(f) = -i[/(4) - 2 /(4-i) + /(*«-»)]• 

(2.106) 

Hence if f(x) is twice differentiable in [xJliV;' 1 ] we get that 


d k } {f) = -UhkffiO ^ some ^ € (x-lj,^ -1 ) 

(2.10c) 


and consequently the scale coefficients in a region of smoothness of f{x) tend to zero as 
Q((/ lfc )2) f or k _► oo. However at a jump discontinuity d k (f) is proportional for the size of 

the jump and thus remains 0(1) as k — *• oo. 

We can obtain data compression by setting to zero all scale coefficients which fall below 
a prescribed tolerance. Let us denote 


f 0 if |d}| < 
d) = 

[ d J if |dj | > £fc, 

(2.11a) 

and 

v L = M _1 • {d L ,...,d 1 ,u°}- 

(2.116) 

Based on the analysis in [H3] we get the following bound on 
case of piecewise-linear interpolation: 

the compression error in the 

L 

max \vi -Vi\<^2e k . 

(2.12a) 

Given any £ > 0 we can take 

e k = 2 k ~ L ~ l e 

(2.126) 

and thus ensure by (2.12a) that 


max |vf — vf | < £• 
0 <i<J L ' 

(2.12c) 
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Example 2. 


Let us consider absolutely integrable functions / in [0,1] 



f€F=L'[ 0,1J, 


and let 


c* = {«?}£. , 4 = (*?., ,*?), 

(2.13a) 

where {x k 

[ are the gridpoints of X k in (2.1); observe that 



-vt’eS 

D 

■i 

II 

7 

(2.136) 


Let v k - be discrete values which are obtained by taking the average over the cells 

in C k of some function / € T = Z ! [0, 1] 


* 

v i =: = pi J ck f dx ' l c ^l = h k, * = 1 < k < L. (2.14a) 

It follows from the additivity of the integral and (2.13b) that v k ~' can be obtained from v k 
by the decimation 

v i = + v 2i)i * = 1 ) • • • ) Jk-1- (2.146) 

Let ( U k v k ){x ) denote any function in Z^O, 1] which satisfies 


(V k ■ n k v k )i 



i 1 )••••> Jk 


e.g. we can take to be the piecewise-constant function 


( 2 . 15 ) 


(1l k v k )(x) = v k for x € c k . (2.16) 

In the context of ENO schemes for the numerical solution of conservation laws we refer to 
Kk as reconstruction from cell-averages and to (2.15) as “conservation” (see [HEOC]). 
Given u* -1 we can approximate v k by 




(T> k -K k - l v k - 1 ) i = 



){x)dx , 


i.e. by first approximating / from v k ~ l by K k - X v h - X and then taking cell-averages of this 
approximation over the finer level k. It follows from the conservation property (2.15) that 
the prediction error e k , 


e i = K - (Vk ■ n k . iv*- 1 ),- , i = i, . . . , j k 


(2.17a) 
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satisfies the relation 


_i + 4 = 9 for i = 1,..., Jk-i 


(2.176) 


(2.186) 


Therefore, if we store 

d k = e k j_i for j = 1,..., Jfc-i (2.18a) 

we can recover the prediction error e k by 

f 4-i = 4 , lx 

J (2.186) 

1 4 = -4 

and thus get a one-to-one correspondence between v k and {d k ,v k 1 }. As in the previous 
example this leads to the multiresolution representation 

v L fiii* {d L ,. . . ,d\u 0 } =: vm 

where the direct MR transform v M = M • can be expressed by the algorithm 
’ DO k = L,...,l 

< v? -1 = |(4-i + 4)> * = ( 2 - 19 ) 

4 = 4-1 - (' Kk-iv k ~ 1 ){x)dx , j = 1,..., Jfc -1 

and the inverse MR transform v L = M -1 • % is given by 

' DO k = 

DO i = 1 

(2.20) 

4-. = ICj + 

, 4 = 2vf _1 - 4-1- 


(2.19) 


( 2 . 20 ) 


Observe that the last statement of (2.20) is obtained from (2.14b). 

In [HEOC] we showed that any interpolation method (2.3) gives rise to a reconstruction 
from cell-averages (2.15) by the following “reconstruction via primitive function” technique: 
Given cell-averages v k = V k f in (2.14a) we calculate 


Ft = F(x k ), F(x) = f X f{y)dy 

J 0 


F k = 0, 4 = 6 fe £4 1 <i<Jk, 

j=i 
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and define 


(R t r*)(x) = —/*(*;/’*), (2.2, 6) 

^ Ft) iS " y ‘"‘“P 01 ' 1110 " “ f ‘he values F‘ = {!?}*„ a , the grid points of x „ in 

d , In ‘ he . Pr ' V,OU * ‘ W0 exam P ,es we have shown how to design MR schemes for discrete 
* ” med fr ° m disc «ti Z ation of functions by point values (Example 1) or by 

coverages (Example 2) in the nested sequence of uniform grids of [0,1], In [H2] and (H31 
we have presented a more genera, class of discrete data which can be repLented by a s cal 
decomposition. This class is characterized by the following notion of nested discretion. 

d“a"o„r ^ tHat 4 ° f ' inear ° Perat ° rS is a -Trence of 


(0 

00 


V k :rZSv\ dim V k = J k , 


Vkf = 0 =» ©*_,/ = 0 


(2.22a) 

(2.226) 


Here T is a space of mappings and V* is a linear space of dimension J t 

data" ‘ 1“ r Sh T hOW l ° ° bUin multiresoluti ™ representation of any discrete 
• , j '{'■ W le ' C he sc, ’ l '-' le "'nipos,ti„n corresponds to the levels of resolution which 
are m reduced in (2.22). This is a very general framework which allows for discretizations 
corresponding to unstructured grids in several space dimensions. 

3. General Multiresolution Representation Schemes 


In this section we consider discrete data which 


is associated with a nested sequence of 


discretion . and show how to design schemes for its multiresolution representation. 

operator ^-'7 h‘ * T* Se<,UenCe diSCre ‘ izati °" with a decimation 

P D k which is a linear mapping from V k = V k {?) onto V k ~' = £> fc _, ( jr) 


D k ~ x 


. yk onto yk-\ 


(3.1a) 


This decimati-on operator is defined as follows: For any , in V* lhere is „ least Qn 
such that V k f = v; the decimation of v is V k ^f e l/* -1 , i.e., 


■V ev k , V = V k f, D k k ~ l v = !>*_„/. 


(3.16) 
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It follows from (2.22b) that 
of the particular / . To see 


is well-defined by (3.1b), i.e. Us definition is independent 
that let us take /i and ft in f such that 


Vkh = v = V k f 2 , 


then by (2.22b) 

0 = Vth _ Vkh = ft(A - h) =*• 0 = ft-> (/■ - w = 

which proves our claim. 

Given v L € V L we can evaluate by repeated decimation 

v k ~ x = k = L,...,l, 


Since (3.1b) implies that 


D k k '\V k f) = Pfc-i/ f° r any / G ^ 


(3.3) 


we get for any / € f for which u‘ = ft/, that u‘ = ft/ for all k in (3.2) We would like 
to stress the point that, as in (2.2b) and (2.14b) in the previous examples, tins dec, mat, on 
is done without explicit knowledge of /. 

Since by (2.22a) V* = ft( F), it follows that V t has a right-, nverse (at least one) w 
we denote by Hk- m n 

n k : V k -* JF, V k K k = Ik, [6A) 

where /. denotes the identity operator in V*. Since (ftu‘) 6 T * an approximation to any 
f € yr for which ft/ = v*, we refer to % as a reconstruct, on of ft; see (2.3)-(2.4) and 

(2. 15)- (2. 16) in the previous examples. L 

Next we show that any sequence of corresponding reconstruct, on operators 

defines a MR scheme for discrete data v L in V L . Starting from v ,n (3.2) we can ge 
approximation to v k by 


rv 

V ~ 


V k {K k - 


We denote , , , /<, 

and refer to it as prediction operator. It follows immediately from taking / = in 

(3.3) and using (3.4) that /?_, is a right-inverse of the decimal, on ft 

BiT'ft, = 'I-- 

We observe that the prediction error e k 


(3.56) 


e fc = v k - - 1 = (/* - 


(3.6a) 


8 



satisfies the relation 


D k k ~ 1 e k = D k k ~ l v k - (D k ~ x P k _ 1 )v k ~ 1 = v^ 1 - = 0 

and therefore it is in the null space of the decimation operator 

e k £N{D k k - l ) = {v | v £ V k , D k ~ l v = 0}. (3.66) 

It follows from (3.1a) that 

dimA/pJ" 1 ) = J k ~ Jk-i (3.7a) 

and hence 

Ar(D i ‘- , ) = span{^}*7''‘-'. (3.7*) 

where {fi k } J j ^[ Jk ~ 1 is any basis of N{D k k ~ l ). Therefore the prediction error e k , which is 
described in terms of J k components in V k , can be represented by its ( J k — Jk- i) coordinates 
d k in (3.7b) 

e " = E ' =■ E * dk ' dk = : G ^ k - (3.8<o 

j'=i 

Here G k denotes the operator which assigns to e k <E Af(D k ~ l ) its coordinates d k in the basis 
{(i k } J j k = r i Jk ~ 1 ; observe that E k Gk is the identity operator in Af(D k ~ l ), i.e. 

E k G k e k = e k for any e k € (3.8 b) 

At this point we can show that there is a one-to-one correspondence between v k and 
{d fc ,u fc-1 }: Given v k we evaluate 

r v k ~ l = D k ~ l v k 

i 

k d k = Gkih-Pt.rD^y ; 

given and d k we recover v k by 

+ E k d k = PLxD k k -'v l ’J r E k G k (h-Pl l D k ,r'> k 
= PLiDi-'v* + (Ik - PL,Dt-')’’ t 
= #*. 

As in the previous examples this shows that 

v iA{d L ,...,d\v°}=: v M , (3-9) 
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where the direct MR transform vm = M ■ v l is given by the algorithm 

' DO k = L,..., 1 

< v k ~ 1 = D k k ~'v k (3.10) 

k d k = G k (I k -P k k _ i D k - 1 )v k =:G?v k 
and the inverse MR transform v L = M -1 • vm can be calculated by 

' DO Jb = 1,...,T 

< (3.11) 

k v k = Pt^v k ~ l + E k d k . 


We remark that in multigrid terminology D k_1 is “restriction” and P k _ x is “prolongation.” 
In signal processing D k k ~ l plays the role of “low-pass filter” while G k , which is defined in 
(3.10), plays the role of “high-pass filter.” 

Example 3. Biorthogonal Wavelets. 


In this example we derive the MR schemes which correspond to the bases of biorthogonal 
wavelets in [CDF]. These MR schemes are obtained from nested discretization of functions in 
T? oc (R) by taking weighted-averages on a nested sequence of uniform grids of R, as follows: 



(£>fc/), = f lk J f( x ) w ^ hk ' ^ dx , OO <1 < oo, 

(3.12a) 

where 

w € T 2 (R) is a weight-funcntion 



roo 

/ w(x)dx = 1 

J — OO 

(3.126) 

and 

x k — ih k , — oo < i < oo, h k — 2~ k h 0 . 

(3.12c) 


In order to obtain a nested sequence of discretization we want to choose w{x ) so that 


N 

{T>k-if)i = X>(Z>*/) 2i-t (3.12 d) 

(= 0 

where ot(, t = 0, . . . , N are real numbers. We observe that since 


f(x) = c — constant => {T ) k /), = c Vi, k 
we have to limit the choice of {c^} by 

N 

= (3.13a) 

(=o 
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From (3.12d) and (3.12a) we get that for any / £ Xj^R) 



this shows that m(x) has to satisfy the functional equation 

N 

w ( x ) = 2j2a ( w(2x - £). (3.136) 

^=0 

This equation has been investigated in [Da] and [CDF] and got the name of “dilation equa- 
tion” in [S]. It is shown there that subject to condition (3.13a) the dilation equation (3.13b) 
has a distribution solution which is determined up to a multiplicative constant and a shift, 
and that u>(x) has a support of size N. Furthermore, if 

X^ q 2 t — (3.13c) 

e t 

then ia(x) is also square integrable. 

We make the choice of w(x) unique by imposing (3.12b) and fixing its support in, say, 
[-N, 0]. With this choice of a weight-function, the sequence of discretization in (3.12) is 
nested and its decimation operator is given by 

(-PaTM. =■' ^2^ev 2i -e = (3.14a) 

t m 

In [Da] and [H3] it is shown that 

= ( — 1 )' +1 ot 2 j— i_i , -oo < i < oo (3.146) 

is a basis of N(D k k ~ x ) in (3.6b). 

We reconstruct the discretization V k in (3.12a) as follows: We take a sequence {p ( } of 
compact support which satisfies 


j He Pit — Ht Pn—\ = l 

(3.15) 

l Hit &tPl+2m — t>m,0 


and define 


(ft*t>)(x) = Y,Vi<p ^ h k ')' 

(3.16a) 

where <p(x) is a solution of the dilation equation 


<p{*) = 'EM 2 * ~ £ ) 

(3.166) 


t 
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which is normalized by 


(3.16c) 


J ip(x)w(x)dx = 1. 

It is easy to see that the corresponding prediction operator P£_ t = Vkllk-i is given by 

( P k-l V )* = : ^2fii-2m v m ( 3 - 17 ) 

m 

Daubechies’ orthonormal wavelets are obtained by imposing the additional condition 

fit = 2 at 


(see [H2] and [H3] for more details). 

We remark that the “fundamental solution” for the dilation equation (3.13b) with 

ot = ^t,o (3.18a) 


is 

w(x) = 6(x), (3.186) 

where 6(x) is the Dirac distribution (see [S]). In this case (3.12a) becomes the point value 
discretization (2.2a) of Example 1. However, since 6(x) is not square integrable, point value 
discretization is excluded from the theory of wavelets. For 

a t = 2^'° ^• 1 ) (3.19a) 


we get 


w(x) = 


1 — 1 < x < 0 

0 otherwise 


(3.196) 


which is square integrable, and the discretization (3.12a) becomes the cell average discretiza- 
tion (2.14a) of Example 2. Observe, however, that the theory of wavelets is for the infinite 
domain R, while our formulation is suitable for both the finite (Example 2) and the infinite 
case (Example 3). Furthermore, unlike the theory of wavelets which uses translates and 
dilates of a single function for both discretization and reconstruction and consequently is 
restricted to uniform grids, our framework of nested discretization allows for general ge- 
ometries. In [H3] and [AH1] we extend the multiresolution representation of cell-averaged 
data in Example 2 to unstructured meshes in bounded domains of R m ,m > 1, by using 
agglomeration of cells to generate a nested sequence of discretization. 

In [ADH] we consider the case where ic(x) in (3.12) is the “hat-function” and show how 
to improve data compression by using adaptive prediction techniques near discontinuities 
and distributions. 
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4. Multiresolution Representation of a 2-Dimensional Array 


In this section we consider functions / 


/ : [0, 1] x [0, 1] — ► R 


which are discretized on the tensor-product grid 


(4.1a) 


x k = {(*?,*;)}&. 


( 4 . 1 4 ) 


A '' 3 ( h k yL Jo f ( Xl ' x ^ w ( hk ') w [ \ k ) dx\dx2 (4.1c) 

where {arf} are the one-dimensional gridpoints in (2.1) and w(x) is a weight-function as in 
Example 3; by (3.18) and (3.19) this includes point value and cell-average discretizations. 
Although this case is covered by the general framework in (3.9)-(3.11), it is convenient here 
to represent the two-dimensional array in (4.1) as the N k x N k matrix A k , and to use tensor- 
product extension of the corresponding one-dimensional operators to get a MR scheme for 
the input A L . 

Let us denote the matrix representation of the various one-dimensional operators by 


D k 1 — ♦ (D) Nk _ l xN k 

P k - 1 ♦ ( P )N k xN k ^ 

Gk = G k (h - PLiDt 1 ) — » (G D ) Nk _ lxNk = G(I - PD) 


(4.2) 


Ek — ♦ (E) NkXNk _ 1 . 

These matrix representations are obtained by taking v k and d k in (3.10)-(3.11) to be column- 
vectors, e.g. 

N k 

v k ~ l = ( D k ~'v k ) t =: J2 D ijV j, 1 <i< N k .! (4.3a) 

i= i 

where by (3.14a) 

Dij = ct2i i-j ; (4.36) 

for simplicity we drop the index k. 

Starting with A L we decimate to get 

A k ~ 1 = DA k D\ k = L,..., 1; (4.4) 

here (•)* denotes the transpose. Given A k ~ l we get an approximation to A k by 


A k « PA k ~ l P* 
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and observe that the prediciton error matrix e k 

e k = A k - PA k ~ x P* 

satisfies 

De k D* = 0. 

The dimension of the null space of the decimation operators here is 

Jk ~ Jk-i = {2Nk-i) 2 - (N k -i) 2 = 3 (Nk-t) 2 

and we store the scale- coefficients d k in three Nk-\ x Nk-i matrices Af, A*, A*. 
Using the matrix identity 

I = PDA EG d , G D = G(I - PD), 

we show that if we take A k and A k ~ x from the sequence (4.4) and define 

A* = G D A k (G D Y , A* = DA k (G D Y, A k = G D A k D* 

then A k can be recovered from A k ~ x and the above by 

A k = PA k ~ 1 P* A EA k E * + PA k E * + EA k P*. 

This follows immediately from the identity 

A k = ( PDAEG D )A k (PDAEG D )* = P(DA k D*)P * 

+ E[G D A k {G D )*]E* A P[DA k {G D Y)E * + E{G D A k D m )P*. 

A 

We conclude from (4.7)-(4.8) that A L has a multiresolution representation Am, 

A l U A m -■ ({ Am} m=1 * - • • > { A m } m=1 » A ) > 

where the direct MR transform is given by 
DO k = L,..., 1 

. A k ~ x = DA k D * 

A* = G D A k {G D Y, A k = DA k (G D Y, A k = G D A k D*, 

and the inverse MR transform is 

(DO k = \,...,L 

\ A k = PA k ~ x P* A EA k E* A PA k E* A EA\P\ 


(4.5a) 

(4.56) 

(4.6) 

(4.7) 

(4.8) 

(4.9) 

(4.10) 

(4.11) 
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In Figure 1, which is taken from [HY], we show the result of applying thresholding with 
tolerance e = 10 -7 to A M in (4.9) where 


— 


! 


1 

i-j 

0 


if i ± j 
if i = j 


(4.12) 


Here we take Nl — 512 and use point value discretization which is reconstructed by a 
sixth-order accurate piecewise-polynomial interpolation with a centered stencil; near the 
boundaries we use one-sided stencils. The rate of compression, i.e. the ratio of (N L ) 2 to the 
number of entries in A M which are above the tolerance e = 10~ 7 , is 8.57. This example is 
taken from [BCR] where it is done with MR schemes which use Daubechies’ wavelets; the 

corresponding rate of compression for wavelets with six vanishing moments is reported to be 
7.33 . 

In Figure 1 we display the results by writing k M in (4.9) as the matrix 





• 

• 

• 

• 



■*! 

- 



A° 


and marking the entries that are larger in absolute value than 10 -7 by a black dot. 

5. Multiresolution Algorithm for Matrix- Vector Multiplication. 

In this section we describe an algorithm to reduce the cost of evaluating the matrix- vector 
multiplication c, 

c = ^>, (5.1) 

where A is a N L x N L matrix which can be thought of as the discretization (4.1) of some 
piecewise-smooth function, and b is any vector of N L components; note that we do not make 
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any smoothness assumptions on 6. As an example let us consider the case where (5.1) is 
obtained from discretization of the integral transform 

c(x) = / a(x,y)b{y)dy, ( 5 - 2a ) 

Jo 

here 

Ci«c(xf), Aij&hL-aix^x!-), bj w b(x]), (5.26) 

and we want to evaluate (5.1) to a specified tolerance which is taken to be of the order of 
the discretization error in (5.2). 

The basic idea is that the product c = Ab, which is being set up as c L = A L b L , has a 
meaningful analog 

c k = A k b k , k = L- 1,...,0 (5.3a) 

corresponding to the k-th grid, and that 

c k = Pc k ~ l + correction, (5.36) 

where P is the matrix representation (4.2) of the one-dimensional prediction operator, and 
the “correction” comes from locations in [0, 1] x [0, 1] where a(x,y ) is still not sufficiently 
resolved on the (k - l)-th grid. To obtain (5.3) let us multiply (4.8) by b k to get 

A k b k = PA k ~\P'b k ) + EA k {E m b k ) + PA k {E*b k ) + EA k 3 (P*b k ), 

from which we see that if we define 

b k ~ l = P*b k , k = L,..., 1 (5-4a) 

and denote 

s k ~ 1 = E*b k , Jb = l,...,L ( 5 * 46 ) 

then we can rewrite the above identity as 

c k = P( M - 1 + E{^ k s k ~ l + A^"- 1 ) + ^(A^*- 1 ). (5.4c) 

Using this observation we get the following algorithm for the approximate (up to a specified 
error) evaluation of the matrix-vector multiplication (5.1). 

(A) Preparation: 

(i) Given A , set A L = A and use the direct MR transform (4.10) to obtain the MR 
representation Am (4-9). 
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(ii) Apply thresholding to Am in order to obtain a sparse representation. 
(B) Multiplication: 

Given any vector b 


(i) Set b L = b 


(ii) Calculate directly 


DO k — L , . . . , 1 
< 6*' 1 = P*b k 
s *- x = E*b k 

c° = A°b° 


(this is done on the coarsest grid). 


(iii) 

I DO for k = 1 , . . . , L 

\ c k = P(c k ~ l + A^s*- 1 ) + E(A k s k ~ l + A^*" 1 ) 

Set c = c L . 


(5.5) 


(5.6a) 


(5.66) 


In the case of pointvalue discretization this algorithm turns out to be identical to the 
multilevel matrix multiplication of Brandt and Lubrecht in [BL]. The algorithm which cor- 
responds to Daubechies’ orthonormal wavelets is identical to the “non-standard form” of 
Beylkin, Coifman and Rokhlin in [BCR]. We remark that (5.5)-(5.6) is a slight generaliza- 
tion of the algorithm which was presented in [HY]. 

Based on the analysis in [BCR] we show in [HY] that if the kernel a(x, y ) in (5.2a) satisfies 

|fl<g(g,y)| < | x ’ * = 0,...,r-l (5.7) 

then outside a diagonal band of width 5, the entries of the matrices {A^}^ =1 are not larger 
than 0(B~ T ), where r is the order of accuracy of the reconstruction technique. Since the 
matrices P and E are banded, we get that the complexity of the algorithm in (5.5)-(5.6) is 

0(N l ) . 

In [HY] we used the above matrix-vector multiplication algorithm to apply the matrix 
in (4.12) to a vector b with “randomly generated” components. Using single precision we 
obtained for the case in Figure 1 a relative residual error ||A6 — c||/||6|| which was 7.52 x 10 -6 
in the norm, and 4.41 x 10 -5 in the l <*, norm. 
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6. Multiresolution Form of Numerical Schemes. 

In this section we consider the numerical algorithm 

u n+1 = Av n + g 


(6.1a) 


which describes either an iterative procedure or a numerical scheme for the solution of an 
initial value problem. Taking the MR transform of (6.1a) we get 

0JJ+ 1 = Mv n+l = {MAM~ l )(Mv n ) + Mg =: A s v n M + g M ; (6.16) 

observe that As = MAM -1 , the multiresolution form of the matrix (= operator) A, is 
different from Am, the multiresolution representation (4.9), where A is treated as a two- 
dimensional array and not as an operator. 

From (3.10) we get that Mv L can be expressed by 

d k = G?{D k k+l ■ ■ ■ D l l - 1 ) • v L =: B k L v L , k = 1, . . . , L 

( 6 . 2 ) 

v° = {D\---D L L ~ l )-v L =:B° l v l . 

From (3.11) we get that M~ 1 vm can be expressed by 


M 1 v>m = Ckd k + CqV°, 

(6.3a) 

k = 1 


Ci = (Pt-, -Pk*')Et . k=l,...,L 

(6.36) 

nL pL pi 

^0 — r L- 1 ’ M> ■ 



where 


Using (6.2)-(6.3) to epxress (6.1b) with g = 0, we get 

d k ' n+ 1 = TS*MACf) ■ d e,n + {B k L AC£) • u 0>n , k = 1, . . . , L 
yo.n+i = '£t =l {B° L ACt') • d e ' n + {B° l AC£) ■ v°’ n 
which shows that 


(6.4a) 
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As = 


B\AC\ 



(6.4b) 


Observe that the block B^AC^ , 1 < k,£ < L, is of size Nk - 1 x N{-\. 


In the Appendix we show that C/' = (B£)* where B e L is given by (6.2) for the dual MR 
scheme; hence 

(6.5) BlACf = BlA(B’ L y 

and each column of (£|;A) is d k , the scale coefficients (6.2) of the k-th column of A, while 
each row of A(Bir is (d e )*, where <r is the column- vector of scale coefficients of the data in 
the ^-th row of A which is obtained by the dual MR scheme. 

We remark that when M in (6.1b) corresponds to Daubechies’ orthonormal wavelets, As 
is identical to the standard form which was introduced in [BCR]. In this case the MR scheme 
is dual to itself, i.e. B £ = B\ and therefore A/ -1 = M* . It was shown in [BCR] that if A is a 
discretization of a Calderon- Zygmund operator, then data compression of each block in (6.5) 
results in a “diagonal” band, the width of which is independent of k and t. Consequently 
applying data compression to A$ results in a finger-like structure of non-zero entries, and 
their total number is 0{Nl log 2 Nj,). 

In [AC] and [ACD] the standard form of [BCR] was extended by (6.1b) to point value and 
cell-average discretizations. It is shown there that, inspite of the lack of symmetry, the rate 
of compression compares favorably to that of orthonormal wavelets: It is about the same in 
the periodic case, but it is significantly better in presence of boundaries. 

In Figure 2, which is taken from [AC], we show the nonzero entries of As, the multires- 
olution form of the matrix in (4.12). As in Figure 1 we use point value discretization which 
is reconstructed by piecewise-polynomial interpolation, however here Ni = 64. 
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The present work extends the scope of application of these algorithms to the general class 
of MR schemes in (3.9)-(3.11). 

7. Multiresolution Application of Operators. 

In this section we describe multiresolution algorithms for the solution of integral equations 
and for the numerical solution of initial-value problems. We shall say that a matrix is T- 
Sparse if it becomes sparse under thresholding, and refer to the number of non-zero elements 
as its T-sparseness. In this language we can express the main result of the previous section 
by saying that the multiresolution form A s of Calderon- Zygmund operators is T-sparse and 
that the T-sparseness of As is 0(Nl log 2 Nl). 

7. A. Integral Transforms and Equations. 

In section 5 we described a multilevel algorithm for matrix- vector multiplication (5.1) 
and showed that if it corresponds to a discretization of an integral transform of the form 
(5.2a) with a kernel which satisfies (5.7), then it can be performed in 0(N L ) operations. We 
can also evaluate c = Ab from its multiresolution form by 

6m = Mb , 

< CM = As'bM, (^-l) 

C = M~ x c m . 

A A 

Since the T-sparseness of A s is 0(Nl log 2 Ni) this is also the cost of the product As&m ; in 
addition we also have to apply the direct MR transform to b and its inverse to c M ■ Hence, 
unless we are in a special situation in which cm and/or &at are also T-sparse, it is more 
efficient to calculate integral transforms by the multilevel algorithm (5.5)-(5.6). 

A situation of this type occurs in a matrix-matrix multiplication C = AB where both 
As and Bs are T-sparse 

C s = MCM" 1 = (MAAT 1 )(MHM" 1 ) = A S B S . (7.2a) 

Of particular interest is the case where the T-sparseness of 

M(A) n M~ x = (A s ) n (7.2 b) 

is uniform in n. In this case (A s ) n for n = 2 m can be computed in m steps of squaring and 
thresholding 

(As) 2 * = [(As) 2 *' 1 ] 2 , k = 1, . . . , m, (7.3) 

where each product above is between matrices with T-sparseness of 0{Nt log 2 Nl). 
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The numerical algorithm (6.1a) can be written as 


v n = {A) n v° + )r{Ay- l g (7.4 a) 

j=o 

or in its multiresolution form (6.1b) as 

v n M = ( As) n v° M + (7.46) 

3=0 

Following [BCR], [EOZ] and [ACD] we use the following algorithm for the fast evaluation of 
(7.4a) for n = 2 m : 


(i) Set 


00 


Calculate 


B = MAM- 1 
C = 1 

DO m times 
< C = tr(C + BC;e) 
B = tr(BB;e) 


(7.5) 


(Hi) 


v n = M- 1 {BMv° + CMg)\ 


here tr(A; e) denotes the truncation operation 


f Aij ^ \Aij\ > £ 

[tr(il;c)]<j = < • (7-6) 

[ 0 if |A,j| < e 

Fredholm integral equations of the second kind are usually solved by iterative procedures 
of the form (6.1a). In this case the gain in efficiency which is offered by algorithm (7.5) is 
based both on the compression of the operator A by 

A — ►fr(A s ;e) (7.7) 

and using the fact that the T-sparseness of (As)" remains constant as n increases. 
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7.B. Initial Value Problems. 

Consider the evolution equation 

d t u + C(x,d x )u = f(x) ,xefi t > 0 
u ( x , 0 ) = uq(x) 


with boundary conditions, where £ is a differential operator. An explicit discretization 
typically takes the form (6.1a) 


v n + l = Av n + g 

(7.9a) 

where 


v] &u(xj,t n ), (/j Ri A(/(ij), 

(7.96) 

t n = nAt and {xj} is a grid in f I. 



It is shown in [EOZ] that using the multiresolution algorithm (7.5) to calculate large- 
time solutions of one-dimensional hyperbolic problems to a fixed predetermined accuracy 
can be reduced from the standard 0((Nl ) 2 ) to 0(AT(log 2 Nl) 3 ). For parabolic equations, 
a standard explicit calculation with complexity 0((Nl) 3 ) can be likewise reduced by (7.5) 
to 0(AT(log 2 Nl) 3 ). The multiresolution algorithm (7.5) in [EOZ] is based on Daubechies’ 
orthonormal wavelets. In [ACD] this algorithm is extended to point value and cell-average 
discretizations. 

As an example let us consider the simplest hyperbolic problem 


Ut — f- Ux — 6 


(7.10a) 


and its solution by the Lax-Wendroff scheme 


,n+i _ 


= V? - - »<"-.) + jV.’V, - 2»r + ».-i) == (Av')i 


(7.106) 


where A — A t/h^. 

The matrix A is a tridiagonal matrix and thus has 3Nl non-zero entries. On the 
other hand As, the multiresolution form of the scheme, has a finger-like structure with 
0(Nl log 2 Nl) non-zero entries. In Figure 3, which is taken from [AC], we show the mul- 
tiresolution form of the Lax-Webdroff scheme for Nl = 64. This shows that unlike the 
application to iterative solution of integral equations where (7.7) results in a “compressed” 
representation of the operator, the gain in efficiency in large-time computation of hyperbolic 
problems is only due to the uniform T-sparseness of powers of As, i.e. while (A) n fills up 
linearly in n, the T-sparseness of (As)” remains 0(Nl log 2 Nl) for all n. 

In the following we describe another application of the multiresolution form (6.1b) which 
uses data-compression of the numerical solution v n and the finite-speed of propagation in 
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hyperbolic problems in order to produce a multiresolution analog to adaptive grids. Let 

r,(»") = {(j, fc)| KV)I>£*} ( 711 ») 

denote the domain in the ( j — /?) plane which contains all the significant scale-coefficients 
of v n , and let f" denote the domain which is obtained by enlarging r e (n”) by adding side- 
neighbors of the cells in r ff (t>”) and allowing for a growth of one scale per time-step where 
needed. Due to the finite speed of propagation in hyperbolic problems 

r ff (u n+1 )cf". (7.116) 

Therefore we can set the components of which are not listed in T" to zero, and evaluate 
the product row(Ts) times only for those rows which are listed in T". The computational 
work can be further reduced by taking into account the T-sparseness of v 1 i . 

This technique can be extended to nonlinear problems. In [LT], [MR] and [BMP] it is 
shown how to derive a multiresolution scheme for the numerical solution of Burgers’ equation 

u t + uu x = vu xx , v > 0 (7.12) 

in which the time-evolution is restricted to the significant scale-coefficients of the numerical 
solution. This numerical scheme is obtained by a Galerkin approach in which the PDE is 
projected on a basis of wavelets. We remark that this Galerkin-type scheme is not suitable 
for the “inviscid” Burgers’ equation (i/ = 0 in (7.12)) in the sense that it generates spurious 
oscillations at shocks, and may even become unstable in some cases - thus some form of 
artificial viscosity is needed. 

In [H4] and [H5] we consider a hyperbolic system of conservation laws 

u t + f(u) x = 0 (7.13a) 

and its numerical solution by any scheme in conservation form 

«;+>=„? -A (fi-fj-i) (7-136) 

where 

fj = f( Vj.K+i , • • • » V ]+K ) for some K ^ (7.13c) 

and / is the numerical flux function. Observe that the computational task here is the evalu- 
ation of the numerical flux function (7.13c) at all the gridpoints. Using the multiresolution 
form of the numerical scheme (7.13b) with respect to cell-average discretization we derive an 
algorithm for the time-evolution of the scale-coefficients, and show that data compression of 
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the numerical solution can be translated into reduction of the number of flux calculations in 
(7.13c). 

In Figure 4 we show the results of [H5] for the multiresolution form of the Lax-Wendroff 
scheme which is applied to the periodic initial value problem for the “inviscid” Burgers’ 
equation. 

u t -f uu x = 0 , — 1 < x < 1 ,f > 0 


u(x, 0) = 2 + sm(7rx) 


(7.14) 


with Jl = 256, CFL = 0.8, and tolerance e = 10 -3 . Figure 4 consists of 3 snapshots 
corresponding to n = 25, 150, 400 time-steps. In the upper part of each snapshot we compare 
the solution of the MR scheme (circles) to the solution of Lax-Wendroff scheme on the finest 
grid (continuous line), which is computed independently. In the lower part of each snapshot 
we display T £ (t; n ) (circles) and its corresponding estimate f " (dots). This is done by drawing 
the symbol at (x^T-i, ^) * n the x — k plane for each (j, k) in the set; note that due to a different 
notation in [H5] k = 0 is the finest grid, and k = L = 5 is the coarsest. In the Table we list 
the efficiency (i.e. the ratio between the fine grid calculation of 256 fluxes over the number 
that we actually had to compute) and the difference E m ,m = 1,2, oo in the corresponding 
norm between the solution of the MR scheme and the independent finest-grid calculation. 

In [BH1] we extend this technique to the numerical solution of 


u t + div f(u) = 0 (7.15) 

on Cartesian grids, and in [AH2] we generalize it further to unstructured meshes where the 
coarser levels of resolution are obtained by agglomeration of cells. 


8. Conclusions. 

In this paper we reviewed recent developments in techniques to represent data in terms 
of its local scale components. These techniques enabled us to obtain data compression by 
eliminating scale-coefficients which are sufficiently small. This capability for data compres- 
sion can be used to reduce the cost of many numerical solution algorithms by either applying 
it to the numerical solution operator in order to get an approximate sparse representation, 
or by applying it to the numerical solution itself in order to reduce the number of quantities 
that need to be computed. 
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Appendix: The Dual MR Scheme. 

In this appendix we describe the MR scheme which is dual to the one in (3.10)-(3.1 1). 
In order to better see the duality we rewrite (3.10)-(3.11) as follows: we express the direct 


MR transform 

by 


where 


We observe that 


vm = M • v l 


DO k = L,...,l 
< v k ~ l = D k k ~ l v k 


d k = Gf v k 

Gg ' =G k (h-p£_ 1 D k k -'). 


(A. la) 


(A.lb) 


(Ale) 


E k d k € M(D k k ~ x ) =► (I k - P k k -xD k k ~ l )E k d k ■= E k d k 
and therefore we can rewrite the inverse MR transform 

v L = M~ 1 vm 


(A.2a) 


as 


( A.2b ) 


DO k = 1 

^ v k = P k _ x v k ~ x + E k d k 

where 

E k = {Ik - P k -\D k k ~ l )E k . (A.2c) 

To simplify our presentation we shall use the matrix representation of the various operators 
and define 

D k ~ l =: {PUT, P k -i=-(D k k ~ l )% G? =:(<)*, E p k =: (G? )*• (A3) 

Observe that Z)£ -1 is a J*_i x matrix, P£_ x is a J k x Jjt_i matrix, Gf is a ( J k - Jfc-i) x J k 
matrix, and E k is a J k x («7^ — J k —\ )• With these definitions we obtain the dual MR scheme 
from (A.1)-(A.2) by applying (•) to all the operators, i.e. the direct MR transform of the 
dual scheme 

(A. 4a) 


is given by 


vr, = M ■ v 1 

M 


f DO k = L , . . . , 1 
v k ~ 1 = D k k ~ 1 v k 
{ d k = Gj?v k 


(A. 46) 
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where Gjf is defined in (A. 3). The inverse MR transform of the dual scheme 

= M - 1 ■ 

is given by 

DO k = 


(A. 5a) 


(A. 56) 


v k = P£-iV k 1 + E£d k 

where E£ is defined in (A. 3). Observe that the dual of the dual is the original scheme. 

It follows from the above definitions that for 1 < k < L 

Ct = (PL, ■ ■ ■ Pt + ')E? = [(E^)•(f^ +, )• ■ ■ ■ (/£_,)T = [Gf (Bt, • • • £>£-')]• = (Blr 

and 


thus 


Cf = Pi-1 ■■■Pi= K^Sr • ■ • (Pt-1 )T = [f? • • • D l l -') ]* = (B»)- ; 


m~' = (My . 


In [H3] we also show that M is associated with discretization f>k and reconstruction 7 Zk 
such that (IZk'Dk) ■ J~ — + T is the adjoint of {TZkDk) • 
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Figure 1 . Multiresolution representation of the 2-D array . 



Figure 2. Multiresolution form of the operator (matrix). 
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